argss <- commandArgs(trailingOnly = TRUE)

print(argss)

pos_f <- argss[1]
seq_f <- argss[2]
prefix <- argss[3]

fill=paste(prefix, 'dup', 'pdf', sep='.')

pos_d=read.table(pos_f, header=T)
seq_d=read.table(seq_f, header=T)


pdf(fill, width=9, height=6)
# par(ann=F, xaxs='i', yaxs='i')

cols=c('blue', 'gray')
# pdf('2.DupRate_plot.pdf')
par(mar=c(5,5,5,10),las=1)
# seq_occ=c(1,2,3)
# seq_uniqRead=c(2288492,693,6)
# pos_occ=c(1,2,3,4,5)
# pos_uniqRead=c(2285791,1845,104,22,3)

# plot(pos_occ,log10(pos_uniqRead),ylab='Number of Reads (log10)',xlab='Frequency',pch=4,cex=0.8,col=cols[1],xlim=c(1,20),yaxt='n')

offse=.2
sub_max=4

seqxmax=max(seq_d$Occurrence)
posxmax=max(pos_d$Occurrence)
seqymax=max(seq_d$UniqReadNumber)
posymax=max(pos_d$UniqReadNumber)
sumseq=sum(seq_d$UniqReadNumber)
sumpos=sum(pos_d$UniqReadNumber)


# xlim=c(0, posxmax+1)
xlimpos=c(0, posxmax+1)
xlimseq=c(0, seqxmax+1)


ylimpos=ceiling(log10(c(.5, posymax)))
ylimseq=ceiling(log10(c(1, seqymax)))

# ylim=ylimpos


# plot(pos_d$Occurrence - offse, log10(pos_d$UniqReadNumber), type='h', ylab='Number of Reads (log10)',xlab='Frequency', col=cols[1],xlim=xlim, ylim=ylim, yaxt='n')

plot(pos_d$Occurrence, log10(pos_d$UniqReadNumber), type='h', ylab='Number of Reads (log10)',xlab='Ruplication Number', main='Read Duplication Rate', col=cols[1], lwd=2, xlim=xlimpos, ylim=ylimpos)

# for (i in 1:nrow(seq_d)) {

# lines(c(seq_d$Occurrence[i] + offse, seq_d$Occurrence[i] + offse), c(0, log10(seq_d$UniqReadNumber[i])), col=cols[2], lwd=2)

# }

# points(seq_occ,log10(seq_uniqRead),pch=20,cex=0.8,col='red')
# ym=floor(max(log10(pos_uniqRead)))
# legend('topright', inset=.05, legend=c('Sequence-base','Mapping-base'), text.col=c(cols[1],cols[2]))
# axis(side=2,at=ylim,labels=ylim)

# labseq=sprintf("%.3f", 10^(seq(ylimseq[2]-1)) / sumseq)
labpos=sprintf("%.2f", pos_d$UniqReadNumber[1:sub_max] / sumpos)

par(cex.axis=.5)
# axis(side=4,at=seq(ylimseq[2]-1), labels=labseq, col=cols[1])
axis(side=4,at=log10(pos_d$UniqReadNumber[1:sub_max]), labels=labpos)


mtext(text = "Reads %", side=4, line = 3)


par(cex.axis=1)
plot(seq_d$Occurrence, log10(seq_d$UniqReadNumber), type='h', ylab='Number of Reads (log10)',xlab='Ruplication Number', main='Read Duplication Rate', col=cols[1], lwd=2, xlim=xlimseq, ylim=ylimseq)

labseq=sprintf("%.2f", seq_d$UniqReadNumber[1:sub_max] / sumseq)
par(cex.axis=.5)
axis(side=4,at=log10(seq_d$UniqReadNumber[1:sub_max]), labels=labseq)

mtext(text = "Reads %", side=4, line = 4)


dev.off()

# s_name=paste(prefix, 'dup', 'rdate', sep='.')
# s_name=paste('.', s_name, sep='')
# save.image(file = s_name)

save.image()
